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ABSTRACT 


Keck-telescope spectrophotometry of the companion of PSR J1810+1744 shows a flat, but asym- 
metric light-curve maximum and a deep, narrow minimum. The maximum indicates strong gravity 
darkening near the Lı point, along with a heated pole and surface winds. The minimum indicates a 
low underlying temperature and substantial limb darkening. The gravity darkening is a consequence 
of extreme pulsar heating and the near-filling of the Roche lobe. Light-curve modeling gives a binary 
inclination 7 = 65.7° + 0.4°. With the Keck-measured radial-velocity amplitude K. = 462.3 + 2.2 
kms~', this gives an accurate neutron star mass Mys = 2.13 + 0.04 Mo, with important implications 
for the dense-matter equation of state. A classic direct-heating model, ignoring the Lı gravitational 
darkening, would predict an unphysical Myns > 3 Mo. A few other “spider” pulsar binaries have similar 
large heating and fill factor; thus, they should be checked for such effects. 


Keywords: pulsars: general — pulsars: individual (PSR J1810+1744) 


1. INTRODUCTION 


PSR J1810+1744 (hereafter J1810) was discovered 
by Hessels et al. (2011) as a P, = 1.7ms, È = 4 x 
1034 I45 ergs~! pulsar in a 350 MHz Green Bank Tele- 
scope (GBT) search of bright unidentified Fermi y-ray 
sources. The dispersion measure DM = 39.7cm~? pc 
indicates a distance of 2.4 kpc (Yao et al. 2017) or 2.0 kpe 
(Cordes & Lazio 2002). It has a double-peaked y-ray 
pulse with flux of 2.3 x 1071! ergcm~?s~+ (Abdo et al. 
2013). J1810 is in a P, = 3.6 hr binary, with a relatively 
large xı = asini = 0.0953 lt-s giving a companion mass 
function 4.2 x 1075 Mo. Thus, for neutron star (NS) 
mass 1.5 Mo, this gives me ~ 0.05 Mo /sini, relatively 
heavy for a black-widow (BW) companion. 

The companion was first identified optically by Breton 
et al. (2013) via 2010 Gemini GMOS-N g-band and i- 
band imaging. Although the light curves were sparse, a 
fit to the data estimated an inclination of i ~ 48° + 7°. 
Schroeder & Halpern (2014) presented more extensive 
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BVR photometry of the bright phases of J1810 using 
the MDM Observatory. They find some disagreement 
with the g photometry of Breton et al. (2013), and quote 
an orbital inclination of i ~ 56° + 3° and a small phase 
shift of the optical maximum from pulsar inferior con- 


junction. 

Here we describe Keck optical spectroscopy and pho- 
tometry of J1810 that allows improved modeling of the 
companion heating pattern. This is important since, in 
addition to direct pulsar y-rays, the companion heating 
can be affected by photons from the system’s intrabinary 
shock (IBS; Romani & Sanchez 2016), by IBS particles 
ducting along magnetic field lines to companion poles 
(Sanchez & Romani 2017), and by heat transfer from a 
global wind (Kandel & Romani 2020, hereafter KR20) or 
general surface diffusion (Voisin et al. 2020). These ef- 
fects distort the light curve, affecting the measurement 
of the system viewing angle and shifting the compan- 
ion center of light (CoL) and radial-velocity center from 
the center of mass (CoM), and can vary from epoch to 
epoch (Kandel et al. 2020), affecting our estimate of the 
NS mass. 


2. OBSERVATIONS 
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Figure 1. Calibrated ugriz light curves from J1810 Keck LRIS spectroscopy (2017), Keck photometry (2020), and GMOS 
photometry (2010). Two cycles are shown and these data are available online as DbF. The model in the first cycle is the best-fit 
direct heating model; the second cycle shows the best-fit HS+GD-+wind model. Residuals to these two models are shown for 
g and 7 in the lower two panels. A few outlier points (marked with black x symbols) are excluded from the fit (see text). See 
Table 1 for the model fit parameters. 


We collected J1810 spectra with the Keck-I 10 m tele- ering 1.1 binary orbits. We used the 5600 A. dichroic, 
scope and LRIS (Oke et al. 1995) for 23 x 600s on the 6001/4000 A blue grism, and the 4001/8500 A red 
UT August 26, 2017 (MJD 57991.285-57991.452), cov- erating, covering ~ 3300-10,500A with dispersions of 
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0.63 A pixel} (blue side) and 1.2 A pixel! (red side). 
The atmospheric dispersion corrector (ADC) allowed us 
to rotate the 1’-wide slit away from the parallactic angle 
(Filippenko 1982) to simultaneously monitor a nearby 
brighter F8 star with known Pan-STARRS2 (PS2) mag- 
nitudes, to monitor the system throughput and wave- 
length solution between frames. In particular, since this 
PS2 monitor star has known and stable magnitudes, we 
integrate the spectra over the SDSS standard ugriz filter 
bands using the sbands IRAF script, and calibrate with 
the catalog magnitudes (converted to the SDSS system 
using the prescription of Finkbeiner et al. 2016) to ob- 
tain light curves with absolute fluxes, up to a possible 
small grey shift from differential slit losses. 

This procedure worked well; the sbands magnitudes 
near maximum brightness are stable and, despite sub- 
stantially decreased throughput toward the end of the 
observation, the fluxes match well in the overlap phases 
(dp = 0.45-0.6). Photometric errors have been esti- 
mated by scaling to the raw spectrum signal. Toward 
minimum brightness, however, the errors are large and 
the cadence too slow for a light curve of suitable quality. 

Thus, we also used LRIS for direct dual-band photom- 
etry, collecting a sequence of g/I images on UT August 
21, 2020 (MJD 59083.38-59083.46) as well as g/R, U/R, 
and g/I images on UT Sep. 17, 2020 (MJD 59100.21- 
59100.38). These had variable 0.8-1.1” seeing. We per- 
formed forced photometry of the companion and a grid 
of nearby PS2 catalog stars in 1.35” diameter apertures. 
The catalog magnitudes, converted to the SDSS scale, 
were used to calibrate the companion light curve and 
estimate systematic errors. In the 0.4” full width at 
half-maximum intensity (FWHM) GMOS-N minimum 
images of Breton et al. (2013), a faint star lies 0.83” 
from the pulsar. We measured fluxes of this star and 
the pulsar companion in all the GMOS-N g/i frames, us- 
ing comparison stars from the PS2 catalog, and subtract 
the estimated contaminator flux in the LRIS photome- 
try apertures. We account for the LRIS/GMOS filter 
and monitor star differences with small shifts (—0.2 mag 
in g, —0.15 mag in i), and find that the points agree well 
with our LRIS-derived fluxes except in the phase range 
0.3 < @g < 0.4, where the GMOS fluxes are ~ 0.3 mag 
brighter than in the LRIS data. This is apparently the 
phase range of a heated spot (see below), and the heat- 
ing may have been stronger in 2010. We also bring 
the LRIS spectrophotometry to the imaging photom- 
etry scale with a gray shift of ~ +0.2 mag, suggesting 
that the monitor-star slit losses were ~ 20% larger than 
those for the companion. 

For the light-curve fitting, we used points with total 
errors (photometric and systematic calibration uncer- 


tainties, added in quadrature) am < 0.15 mag. However, 
we dropped two g points in twilight and one i point taken 
as the source was setting, and excluded the sbands pho- 
tometry at 0.1 < pg < 0.4 (low signal-to-noise ratio and 
large Ad). Altogether we have 23 u, 60 g, 34 r, 51 i, 
and 17 z measurements (Fig.1). Both our new light 
curve and the MDM BVR photometry of Schroeder & 
Halpern (2014) show a rather broad, flat maximum, with 
a slight brightening to later phases and little color vari- 
ation. Our Keck imaging with its fine cadence shows a 
narrow light-curve minimum, reaching g > 25mag. The 
night-phase light curve is significantly asymmetric, with 
the minimum at ¢p < 0.25. 

The extinction in this direction is estimated from the 
three-dimensional dust maps of Green et al. (2018), 
reaching its maximal E(g — r) = 0.12 + 0.02 mag (Ay = 
0.39 + 0.07 mag) by 1.1 kpc. 

Balmer absorption lines dominate the companion 
spectrum through the bright half of the orbit (spectral 
class ~A2 at maximum brightness). Accordingly, we 
initially measure the radial-velocity (RV) amplitude by 
cross-correlation with an A2 template. The correlation 
during the bright “day” phases gives RV uncertainties 
as small as 5kms~!. Although the effective temper- 
ature drops dramatically, the correlation persists with 
plausible velocities (with large uncertainty) for almost 
all night-phase spectra. We have chosen to focus on 
the measurements with the best correlation coefficient 
(R > 10). This gives 15 RV measurements covering 
phase ¢g = 0.45-1.05. We have removed the nominal 
systemic velocity. A simple sinusoidal fit to these data 
gives Kobs = 426.9+3.4kms~! with x?/DoF = 2.3 (Fig. 
2). 


3. PHOTOMETRIC/RADIAL-VELOCITY FITTING 


Our fits are performed with an outgrowth of the 
ICARUS light-curve model (Breton et al. 2012) with ad- 
ditions described by KR20 and Kandel et al. (2020). 
An additional update replaced the simplified limb- 
darkening laws in the base code with the more detailed 
limb-darkening coefficients computed by Claret & Bloe- 
men (2011) for two models, the ATLAS and PHOENIX 
atmospheres. We generally find that the ATLAS coeffi- 
cients perform better. Note that gravity and limb dark- 
ening serve to rescale the local fluxes; we take care to 
integrate the emergent flux to determine the total non- 
thermal pulsar heating Ly and the thermal base emis- 
sion (characterized by the night-side temperature Ty). 

All fits find that the companion is very close to Roche- 
lobe filling, so we set the fill factor at fı = 0.99. We 
also find that, when left free, the extinction prefers val- 
ues larger than estimated from the dust maps. We have 
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Figure 2. Radial-velocity (RV) measurements for J1810. 
The first cycle shows 15 high-correlation (R > 10) measure- 
ments against an A2 template (velocities available on-line 
as DbF). These are fit to a simple sinusoid (dashed curve) 
with the fit residuals indicated at bottom. The curve at 
bottom shows the sinusoid residuals expected for the best 
HS+GD+wind photometric fit; the fit residuals follow this 
curve. The second cycle includes low-significance velocity 
correlations from the night half of the orbit. The solid line 
is again the best-fit photometric model. For this model, 
the phase spectral templates provide 17 measurements with 
R > 10, decreased residuals (shown in the bottom section), 
and an improved x?/DoF = 1.2 fit. 


chosen to fix this at Ay = 0.5mag for our basic fits, 
about 1.5g larger than the dust estimate, and note be- 
low when a free fit gives a significantly different value. A 
standard direct heating (DH) fit gives i ~ 53°; similar 
results were obtained by Schroeder & Halpern (2014). 
The cause is the model’s attempt to produce a flat max- 
imum via low inclination, 7. With our newly measured 
K, this value of i implies an unphysical My > 3 Mo 
(Table 1). 

Asymmetry of the light-curve maximum and bluer col- 
ors leading the peak indicate that there is extra heat 
at this phase. We model this with a Gaussian hot 
spot (HS) of radius ry, centered at (Ons, ns), heated 
o (1+ Ans) times the local star temperature. In the 
picture of Sanchez & Romani (2017), such hot spots 
are caused by relativistic particles ducting from the pul- 
sar wind and IBS to companion magnetic poles. Note 
that the excluded GMOS data points near the HS phase 


(dp ~ 0.4) indicate that this phase was brighter and 
bluer in 2010; we speculate that the particle flux duct- 
ing to the HS varies between epochs (as seen for PSR 
J2339—0533 by Kandel et al. 2020), and was stronger at 
the GMOS epoch. 

We next amend the gravity darkening (GD; Espinosa 
Lara & Rieutord 2012). While the standard ICARUS 
code applies GD to the underlying companion, it does 
not apply to the heated companion face. J1810 has a 
Roche lobe fill factor fı ~ 1 (giving very low g near 
the nose) and temperatures on the heated face exceed 
10,000 K, so GD effects can be substantial, We use a 
simple prescription T” = T(g/go)°. For the effectively 
radiative photosphere of the strongly heated day face 
we use 3 = 0.25 (as usual 6 ~ 0.08 applies for the low- 
T convective atmosphere on the back side), with the 
scaling go taken from the dawn equator point. The pre- 
GD temperature distribution includes such effects as hot 
spots and surface heat redistribution, if any. 

To produce the broadening and slight gradient across 
the maximum, some flux must be moved to the trailing 
side. We find that a simple global wind, with heat ad- 
vected along latitude lines (see KR20 for model details), 
greatly helps in matching the maximum shape. An al- 
ternative is to invoke heat diffusion away from the com- 
panion nose, as described by Voisin et al. (2020). This 
does broaden the peak, but without the asymmetry re- 
produced with the HS+GD-+wind model. Table 1 shows 
the nested light curve fits as these effects are turned on; 
the large x? decreases indicate that the added parame- 
ters are highly significant. 

With these effects included, we have a much better 
companion light-curve fit; wind and GD flatten the day- 
face temperature, while a hot spot reproduces the ob- 
served peak asymmetry. The very narrow deep min- 
imum is not fully captured. The model minimum is 
sensitive to the limb-darkening coefficients: use of the 
PHOENIX model-derived coefficients produces signifi- 
cantly worse fits than those with the ATLAS model coef- 
ficients, possibly due to a different treatment at the low 
temperatures found on the ‘night’ side. Although not 
employed here, a further small reduction in x? can be 
found by decreasing the limb darkening at low temper- 
ature. Thus, additional physical ingredients are needed 
and might be probed with more data and improved at- 
mosphere modeling. 

We next compare the Keck spectroscopy with the 
photometric-fit models. As emphasized by Linares et al. 
(2018) and discussed by KR20, different species’ temper- 
ature sensitivities cause varying line equivalent widths 
(EWs) across the face of the companion. As the model 
parameters change, the surface heating changes and the 
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Table 1. Light Curve and Template Spectral Fit Results for J1810 


Parameters DH HS HS+GD HS+GD+wind HS+GD/Ay * HS+GD-diff 2HS+GD 
i (deg) 52.9108 56.619. 61.7497 wee 66.1705 68.0419 62.7404 
Ly /10°4 (erg/s)  7.04+9:23 6.315517 6.1970 33 6.007906 6.757034 10.970-23 5.9140-07 
Tw (K) 33201508 3440780 3550720 3470 + 25 361073? 24301132 3560732 
dkpe 3.21+ 0.07 3.10+0.06 3.02 + 0.03 3.03 + 0.01 2.91 + 0.02 3.06 +0.05 3.06 + 0.02 
Ons (deg) z 108.1133 112.6421 100.6413 118.4729 117.5735 106.6132 
dhs (deg) - 43.3776 46.0788 2277s 51.4439 53.8736 30.6738 
Ahs - 1.9708 15203 1.2+0.1 1.3+0.1 21+0.4 1.2+ 0.1 
Ans2 s - - - - 0.80.1 
Ths (deg) - 19:7434 18.8736 257542 18.5778 16.6137 20.2713 
€ < - —0.104 + 0.005 - - - 
diff = = = - = 8.9+1:8 = 
x2/DoF 3465/179 1750/175 886/175 248/174 535/173 523/174 270/174 
Kcom (km/s) 474.9+2.0 4730419 468.3 + 2.0 462.3 + 2.2 462.5 + 2.6 463.2+2.1 465.04 2.1 
Mys (Mo) 3.42+0.09 2.95 +0.07 2.45 + 0.05 2.13 + 0.04 2.11 + 0.05 2.03 +0.08 2.33 + 0.04 
Mc (Mo) 0.101 + 0.002 0.087+ 0.002 0.073+0.001 0.065 0.001 | 0.064+ 0.001 0.062 0.003 0.070 + 0.001 
x? 16.53 16.08 16.69 18.28 19.41 18.38 16.52 
x2/15 DoF 1.10 1.07 1.11 1.22 1.29 1.23 1.10 


* Fit 8 = 0.46 + 0.02, Ay = 0.60 + 0.03 mag 


line EWs vary. Since our spectra are strongly Balmer 
dominated, except at @g ~% 0.25, a simple temperature- 
dependant weight for the Balmer EW allows a fit to the 
A2 cross-correlation measurements of the RV parame- 
ters, while allowing the heating model to vary. We do 
not fit RVs simultaneously with the photometric data; 
however, the RV fits for the CoM velocity and the corre- 
sponding masses are marginalized over the parameters of 
every tenth model from the last 1000 of the photometric 
nested sampling chain (Buchner et al. 2014), sampling 
2o uncertainties. Thus, we obtain errors on the CoM K 
(ox) and total uncertainties on the component masses, 
including the uncertainties in the photometrically deter- 
mined parameters (e.g., inclination 7). The photometric 
parameter contribution @pnot can be isolated by propa- 
gating the contribution of ox and subtracting in quadra- 
ture: Ophot / MNs = (0 Mys /Mys}? — (30xK/K)?]!/?. 

Of course, the companion presents multi-temperature 
spectra, varying with phase, so the initial A2 templates 
are imperfect representations. However, once a photo- 
metric model has fit parameters, we can use this model 
to compute in detail, for the model CoM K, the ex- 
pected spectrum collected during each observational in- 
tegration, including all surface temperature and log g 
effects and the Doppler shifts associated with each re- 
gion of the surface. These integrated spectral models, 
computed with ~ 10 times the data spectral resolution, 
should better represent the line strengths and Doppler 
distortions of the companion spectra than the original 


A2 templates. We can thus cross-correlate the observed 
spectra with these multi-temperature templates to ob- 
tain corrections to the model’s CoM RV at each phase. 
Indeed, with these templates 17 spectra now have cor- 
relation R > 10. The cross-correlation residuals to the 
model RVs can be fit with a simple sinusoid to give im- 
proved measurements of the CoM RV K, with smaller x? 
and ox uncertainties. Adding in quadrature the opnot 
contributions measured with the A2 template fits, we 
obtain the final mass error estimates. These are re- 
ported in the bottom section of Table 1, and the velocity 
residuals for spectral template fitting of the best-fit pho- 
tometric model are shown in the lower panel of Figure 
2 (second cycle). 

Note that as the heating model varies, the spectral 
template and hence cross-correlation RV for each ob- 
servation will change. Recomputing the model spectra 
and remeasuring the velocities at each fit step would be 
computationally prohibitive. Thus, our hybrid scheme, 
fitting for a fixed A2 template to capture the RV cor- 
relations with the (photometrically determined) heating 
parameters, while using the composite template model 
spectra for the final CoM RV amplitude Kcom, allows 
us to capture both heating model systematics and an 
accurate spectral representation for the final velocities. 
The RV changes from the A2-determined estimates are 
only a few kms~', so this partition of the fits should be 
robust. 
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As noted, with the DH model assumed in most other 
spider binary studies, the large observed velocity am- 
plitude of J1810 gives an unphysically high NS mass. 
Successively including companion HS, GD, and wind ef- 
fects moves the CoL back toward the CoM, and (as seen 
above) increases the best-fit inclination. This substan- 
tially decreases the inferred NS mass. The y?/DoF of 
the RV fit with the HS+GD-+wind fit is also lower than 
that of a simple sinusoidal RV fit (Fig. 2), implying that 
the heating distortions of the model are reflected in both 
the photometric and spectroscopic data. However, the 
x? differences between the various heating models are 
not significant in the RV fits; all models are acceptable 
and the HS+GD-+wind model has an RV-fit probabil- 
ity only 27% lower than that of the DH model. Model 
selection is thus based on the photometric fit. 


4. MODEL-FIT COMPARISONS 


We have demonstrated that a simple direct (photon) 
heating model does not adequately describe the light 
curve of J1810. However, with the effects of a localized 
hot spot, global winds, and an improved treatment of 
gravity and limb darkening, we get a reasonable photo- 
metric fit with y?/DoF = 1.4 dominated by local un- 
modeled effects at minimum brightness. We have at- 
tempted to see if other models can do as well, but the 
(HS+GD-+wind) model remains the best fit. Before de- 
scribing these models, it is interesting to note the com- 
mon features of all viable fits. First, the companion is 
very close to filling its Roche lobe. Next, all successful 
models require that the companion “nose” near L1 be 
appreciably gravity darkened. 

All successful models also require extra heating in a 
spot past the dawn terminator, giving excess flux at 
phase ¢g ~ 0.4. Such features suggest that the com- 
panion supports a large global field to channel parti- 
cles to magnetic poles. Such spots are commonly found 
on the redback systems with larger (> 0.08 Mo), core- 
fusion-supporting companion masses, but are not usu- 
ally prominent in true BWs. It is interesting that 
J1810’s companion is the most massive known among 
the BWs (although apparently below the core-fusion 
threshold). We have weak evidence that the heating flux 
varies, with larger values during the 2010 GMOS obser- 
vations. It will be interesting to determine whether, as 
for some redbacks (see Kandel et al. 2020), the heated 
spots vary in both flux and position. Precision multi- 
color light curves following such variation could further 
refine the underlying heating pattern. 

Another feature common to all the fits is a large ra- 
diation (DH) flux, Ly > 6 x 10% ergs~?. At first sight, 
this might seem in conflict with the observed Fermi flux, 


which for an isotropic emitter at our fit distance implies 
L., = 2.5 x 10° ergs~ +. However, modern high-altitude 
y-ray emission models, from the outer magnetosphere 
or near wind-zone, direct the radiation toward the spin 
equator (and companion) more than the isotropic Ly 
estimate in Table 1. Thus, the true Ly is lower by 
a model-dependant beaming factor (see Draghis et al. 
2019). Also, one might worry that Ly exceeds the 
spin-down power È = 4 x 10°44; ergs~!. However, 
the beaming correction helps, and one should remem- 
ber that with the large mass inferred here we expect 
I45 ~ 2-2.5 for the stiff equations of state that allow 
2 Mo NSs. 

Table 1 shows the dramatic fit improvements as HS, 
GD, and wind effects are added, justifying our rela- 
tively complex model. However, one can imagine other 
physical effects producing a similar heat distribution; 
we show three such models in the second section of 
the Table, with similar parameter count to our best-fit 
(HS+GD-+wind) model. All require an HS at nearly the 
same location. For example, if we free gravity darkening 
and extinction values (HS+GD/Aj,), Av increases fur- 
ther above the dust estimate, the fit 6 = 0.46 nears the 
largest plausible values (Claret & Bloemen 2011), and 
yet x? is still two times that of our best-fit model. A sim- 
ilar x? can be achieved by replacing the global wind by 
heat diffusion away from local maxima (HS+GD-+diff). 
While this diffusion broadens the peak, it does not du- 
plicate the peak gradient. Notice that other parameters 
(e.g., inclination 7, HS location) are very similar for these 
models. With a similar heating pattern, it is unsurpris- 
ing that best-fit pulsar masses are within lo. 

The last alternate model posits that the companion 
magnetic field is dipolar, with an antipodal HS in the 
southern hemisphere, which serves to add flux, broad- 
ening the light-curve maximum. This spot is actually 
on the companion day side, below the equator, and the 
best-fit model drives to slightly lower inclination i to de- 
crease its contribution to the light curve. Although this 
has the best x? of the alternative models with similar pa- 
rameter count, we do not favor it; the second (day-side) 
spot actually has a lower fit flux (we expect the parti- 
cle ducting there to be stronger) and the suppressed i 
increases the fit mass by 3.40. 


5. MASS IMPLICATIONS AND CONCLUSIONS 


Including the improved RV measurements with com- 
posite template spectra, the statistical error on the NS 
mass of our best-fit model, including all spectroscopic 
and light-curve effects, is quite small at 0.04 Mo. For 
related single HS physical models, the mass changes are 
small (~ 1c) but the x? values are substantially worse. 
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The 2HS model gives a higher though possibly accept- 
able y?, but requires a 0.2. Mọ mass increase. We con- 
clude that J1810 is a heavy NS at 2.13+0.04 Mo, with 
a robust lower limit; poorer, but possibly acceptable, 
models allow masses as large as 2.3 Mo. 

As noted above, J1810 has the largest companion mass 
among the true BWs. Whether this represents an early 
stage of companion evaporation or an unusual early ter- 
mination of the mass-transfer phase is unclear. The high 
heating luminosity and large fill factor raise the possibil- 
ity that the pulsar may transition back to an accretion 
phase. If so, it must be on a decades timescale as J1810 
does not show strong variability in the Fermi flux record. 

The existence of heavy (~ 2 Mọ) NSs was assured 
by the Shapiro delay measurement of PSR J1614—2230 
(Demorest et al. 2010), although subsequent observa- 
tions have refined and lowered the mass to 1.908 + 
0.016 Me (Arzoumanian et al. 2018). Two other mass 
measurements have proved very influential in the liter- 
ature: PSR J0348+0432 at 2.01 + 0.04 Mọ (Antoniadis 
et al. 2013) and the recent Shapiro-delay measurement 
of PSR J0740+6620 at 2.144959 Mo (Cromartie et al. 
2020), with all uncertainties 1o. By making a good- 
precision measurement of Myg for J1810, we improve 
the lower bounds on Mmax. Indeed, this is the first in- 
dividual object for which a 30 lower bound on the mass 
exceeds 2 Mo. We can supplement these masses with 
the measurements for two other spider binaries from 
KR20: PSR J1959+2048 at 2.18 + 0.09 Mọ and PSR 
J2215+5135 at 2.24 + 0.09 Mo (although that analysis 
did not include GD effects so the actual inclinations and 
masses could be somewhat lower). Figure 3 shows the 
mass uncertainty ranges for these objects. 

Several approaches can be used to estimate Max. 
One option is to model the full distribution of (binary) 
NS masses and see if an upper cutoff is required; Als- 
ing et al. (2018), for example, determine that Mmax 
is within a lo range of 2.0-2.2 Mọ. Here we only at- 
tempt to determine a lower bound to Max, so with 
individual source probability density functions (PDFs) 
P; = (2ro?)—1/2¢e-0-5[(m—m)/o:]” we can form the joint 
probability of getting the set of measurements {d;} = 
{mi,oi} when M < Mmax as I; fo" P;dm. This is 
also the Bayesian probability P(Mmax|d;) for a flat prior 
with a hard cutoff, O(m — Mmax). These are shown in 
the lower panel of Figure 3. 

As demonstrated in this paper, we have marginal- 
ized over all parameters in determining the spider mass- 
estimate uncertainties. We have also explored alterna- 
tive models; in most cases these (worse-fitting) mod- 
els require higher masses. Since the spider binary mass 
estimates rely on the heating model, we cannot claim 
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Figure 3. Mass estimates for heavy NSs. The dashed curves 
show three radio-selected pulsars with white dwarf (WD) 
companions, measured from pulse timing (supplemented by 
WD atmosphere modeling for J0348). J1614 (amplitude de- 
creased by a factor of 2 for plot) no longer significantly adds 
to the Mmax lower bound. The solid curves show three spi- 
der binary mass estimates, relying on companion spectropho- 
tometry. The well-determined J1810 mass range is shaded. 
The bottom panel shows the cumulative probability distri- 
butions for M < Mmax, for the radio objects, J1810 alone, 
and all six pulsars. Lower bounds (10, 2c, and 30) are listed 
for these distributions. 


that they are free of systematic effects. However, we 
have attempted to be conservative — and as more spi- 
der pulsars are measured with high masses and increas- 
ing accuracy, we should not ignore their contribution 
to constraints on Mmax and thus on the dense-matter 
equation of state. Taken at face value, the uncertainties 
in recent spider measurements are sufficiently small to 
significantly improve the bounds. We can now say with 
high (30 one-sided lower bound) statistical confidence 
that Mmax > 2.12 Mo, and that at ~ lo significance 
Mmax > 2.24Mo is preferred. Additional measure- 
ments of spider binaries, especially the extreme BWs, 
will tighten (and likely slightly raise) this lower bound. 

In summary, we have found a good photomet- 
rict+spectral model for BW J1810. This relies on sub- 
stantial GD on the heated (day side), especially near 
the companion nose. Also important are HS and limb- 
darkening effects, particularly near binary minimum. 
These effects all serve to increase the best-fit inclina- 
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tion į and to lower the inferred NS mass. Other binaries 
should be checked for this GD effect, but few have the 
large fill factor and high T that make GD so strong in 
J1810. However, even after modeling this effect, J1810’s 
NS mass is large, and alternate (worse-fitting) models 
tend to be even heavier. With our excellent fit preci- 
sion, J1810 provides for the first time a lower limit on 
an NS mass that is greater than 2 Mọ at > 30 con- 
fidence. This seems robust to any residual systemat- 
ics and should thus be important for discussions of the 
dense-matter equation of state. 
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